﻿#pragma once

#include <learn/geom/bounding_box.h>
#include <learn/geom/geom_bspline_curve.h>

namespace xi
{

namespace alg
{

using namespace xi::math;
using namespace xi::geom;

template <std::size_t N> struct Alg_Parameter
{
    using Vec = VecT<N, double>; // P
    using VecP = VecT<N - 1, double>;

    /**
     * @brief 均匀参数化 uk = k / n, k = 0, 1,..., n
     *
     * @param[in] points
     * @return std::vector<double>
     */
    static std::vector<double> equally(const std::vector<Vec> &points)
    {
        std::size_t n = points.size();
        std::vector<double> para(n);
        double step = 1.0 / (n - 1);
        for (std::size_t i = 0; i < n; i++)
            para[i] = step * i;
        para[n - 1] = 1.0;
        return para;
    }

    /**
     * @brief 弦长参数化
     *
     * @param[in] points
     * @return std::vector<double>
     */
    static std::vector<double> chord(const std::vector<Vec> &points)
    {
        std::size_t n = points.size();
        std::vector<double> para(n);
        para[0] = 0;

        double total = 0;
        for (std::size_t i = 1; i < n; i++)
            total += points[i].distance(points[i - 1]);

        for (std::size_t i = 1; i < n - 1; i++)
        {
            double di = points[i].distance(points[i - 1]);
            para[i] = para[i - 1] + di / total;
        }
        para[n - 1] = 1.0;
        return para;
    }

    /**
     * @brief 向心参数化
     *
     * @param[in] points
     * @return std::vector<double>
     */
    static std::vector<double> centripetal(const std::vector<Vec> &points)
    {
        std::size_t n = points.size();
        std::vector<double> para(n);
        para[0] = 0;

        double total = 0;
        for (std::size_t i = 1; i < n; i++)
            total += std::sqrt(points[i].distance(points[i - 1]));

        for (std::size_t i = 1; i < n - 1; i++)
        {
            double di = std::sqrt(points[i].distance(points[i - 1]));
            para[i] = para[i - 1] + di / total;
        }
        para[n - 1] = 1.0;
        return para;
    }

    /**
     * @brief 均匀参数化 uk = k / n, k = 0, 1,..., n
     *
     * @param[in] points
     * @return std::vector<std::vector<Vec2>>
     */
    static std::vector<std::vector<Vec2>> equally(const std::vector<std::vector<Vec>> &points)
    {
        std::vector<std::vector<Vec2>> para(points.size(), std::vector<Vec2>(points[0].size()));

        double ustep = 1.0 / (points.size() - 1);
        double vstep = 1.0 / (points[0].size() - 1);

        for (std::size_t i = 0; i < points.size(); i++)
        {
            double u = ustep * i;
            for (std::size_t j = 0; j < points[0].size(); j++)
            {
                double v = vstep * j;
                para[i][j] = Vec2{u, v};
            }
            para[i].back() = Vec2{u, 1.0};
        }
        para.back().back() = Vec2{1.0, 1.0};
        return para;
    }

    /**
     * @brief 弦长参数化
     *
     * @param[in] points
     * @return std::vector<std::vector<Vec2>>
     */
    static std::vector<std::vector<Vec2>> chord(const std::vector<std::vector<Vec>> &points)
    {
        std::vector<std::vector<Vec2>> para(points.size(), std::vector<Vec2>(points[0].size()));
        for (std::size_t i = 0; i < points.size(); i++)
        {
            double total = 0;
            for (std::size_t j = 1; j < points[0].size(); j++)
                total += points[i][j].distance(points[i][j - 1]);

            para[i][0][1] = 0;
            for (std::size_t j = 1; j < points[0].size(); j++)
            {
                double di = points[i][j].distance(points[i][j - 1]);
                para[i][j][1] = para[i][j - 1][1] + di / total;
            }
            para[i].back()[1] = 1.0;
        }

        for (std::size_t i = 0; i < points[0].size(); i++)
        {
            double total = 0;
            for (std::size_t j = 1; j < points.size(); j++)
                total += points[j][i].distance(points[j - 1][i]);

            para[0][i][0] = 0;
            for (std::size_t j = 1; j < points.size(); j++)
            {
                double di = points[j][i].distance(points[j - 1][i]);
                para[j][i][0] = para[j - 1][i][0] + di / total;
            }
            para.back()[i][0] = 1.0;
        }
        return para;
    }

    /**
     * @brief 向心参数化
     *
     * @param[in] points
     * @return std::vector<std::vector<Vec2>>
     */
    static std::vector<std::vector<Vec2>> centripetal(const std::vector<std::vector<Vec>> &points)
    {
        std::vector<std::vector<Vec2>> para(points.size(), std::vector<Vec2>(points[0].size()));
        for (std::size_t i = 0; i < points.size(); i++)
        {
            double total = 0;
            for (std::size_t j = 1; j < points[0].size(); j++)
                total += std::sqrt(points[i][j].distance(points[i][j - 1]));

            para[i][0][1] = 0;
            for (std::size_t j = 1; j < points[0].size(); j++)
            {
                double di = std::sqrt(points[i][j].distance(points[i][j - 1]));
                para[i][j][1] = para[i][j - 1][1] + di / total;
            }
            para[i].back()[1] = 1.0;
        }

        for (std::size_t i = 0; i < points[0].size(); i++)
        {
            double total = 0;
            for (std::size_t j = 1; j < points.size(); j++)
                total += std::sqrt(points[j][i].distance(points[j - 1][i]));

            para[0][i][0] = 0;
            for (std::size_t j = 1; j < points.size(); j++)
            {
                double di = std::sqrt(points[j][i].distance(points[j - 1][i]));
                para[j][i][0] = para[j - 1][i][0] + di / total;
            }
            para.back()[i][0] = 1.0;
        }
        return para;
    }

    /**
     * @brief OBB 包围盒投影参数化
     *
     * @param[in] points
     * @return std::vector<VecP>
     */
    static std::vector<VecP> projection(const std::vector<Vec> &points)
    {
        auto bbox = BoundingBox<N, double>::OBB(points);

        auto &min = bbox.min;
        auto &max = bbox.max;
        auto &params = bbox.params;
        for (int i = 0; i < points.size(); i++)
        {
            auto &p = params[i];
            for (std::size_t j = 0; j < N - 1; j++)
                p[j] = (p[j] - min[j]) / (max[j] - min[j]);
        }

        return bbox.params;
    }
};

template <std::size_t N> struct Alg_Knot
{
    using Vec = VecT<N, double>; // P

    /**
     * @brief 均匀节点向量，n 个控制点的 p 次样条
     *
     * @param[in] n
     * @param[in] p
     * @return std::vector<double>
     */
    static std::vector<double> equally(std::size_t n, std::size_t p)
    {
        std::vector<double> knots(n + p + 1);
        for (std::size_t i = 1; i <= p + 1; i++)
            knots[i - 1] = 0;
        double step = 1.0 / (n - p);
        for (std::size_t i = p + 2; i <= n; i++)
            knots[i - 1] = step * (i - p - 1);
        for (std::size_t i = n + 1; i <= n + p + 1; i++)
            knots[i - 1] = 1;
        return knots;
    }

    /**
     * @brief 平均节点向量，用于插值算法
     *
     * @param[in] params
     * @param[in] p
     * @return std::vector<double>
     */
    static std::vector<double> average(const std::vector<double> &params, std::size_t p)
    {
        // 注意确保参数值归一化
        std::size_t n = params.size();
        std::vector<double> knots(n + p + 1);
        for (std::size_t i = 1; i <= p + 1; i++)
            knots[i - 1] = 0;

        for (std::size_t i = p + 2; i <= n; i++)
        {
            double s = 0;
            for (std::size_t j = 1; j <= p; j++)
                s += params[i - p - 2 + j];
            knots[i - 1] = s / p;
        }

        for (std::size_t i = n + 1; i <= n + p + 1; i++)
            knots[i - 1] = 1;
        return knots;
    }

    /**
     * @brief 均匀节点向量，n 个控制点的 p 次样条，用于逼近算法
     *
     * @param[in] params
     * @param[in] n
     * @param[in] p
     * @return std::vector<double>
     */
    static std::vector<double> uniform(const std::vector<double> &params, std::size_t n, std::size_t p)
    {
        // 注意确保参数值归一化
        std::size_t m = params.size();
        double d = 1.0 * m / (n - p);
        std::vector<double> knots(n + p + 1);
        for (std::size_t i = 1; i <= p + 1; i++)
            knots[i - 1] = 0;

        for (std::size_t j = 1; j < n - p; j++)
        {
            int i = j * d;
            double alpha = j * d - i;
            knots[j + p] = params[i - 1] * (1 - alpha) + params[i] * alpha;
        }

        for (std::size_t i = n + 1; i <= n + p + 1; i++)
            knots[i - 1] = 1;
        return knots;
    }

    /**
     * @brief 通过节点插入和升阶使两个 B-spline 曲线的节点和次数相同
     *
     * @param[in] c1
     * @param[in] c2
     */
    template <bool Rat = false>
    static void compat(std::shared_ptr<Geom_BSplineCurve<N, Rat>> c1, std::shared_ptr<Geom_BSplineCurve<N, Rat>> c2)
    {
        std::size_t d = std::max(c1->degree(), c2->degree());

        c1->elevate(d);
        c2->elevate(d);

        // 对齐两个曲线的第一个和最后一个节点
        c1->reparam(c2->first_param(), c2->last_param());

        auto knots1 = c1->knots();
        auto knots2 = c2->knots();

        std::vector<double> ins1;
        std::vector<double> ins2;

        std::size_t i = d + 1, j = d + 1;
        while (i < knots1.size() - d - 1 || j < knots2.size() - d - 1)
        {
            auto d = knots1[i] - knots2[j];
            if (std::abs(d) < Numeric::half)
            {
                // 节点重合，跳过
            }
            else if (d < Numeric::half)
            {
                ins2.push_back(knots1[i]);
                knots2.insert(knots2.begin() + j, knots1[i]);
            }
            else if (-d < Numeric::half)
            {
                ins1.push_back(knots2[j]);
                knots1.insert(knots1.begin() + i, knots2[j]);
            }

            if (i < knots1.size() - d - 1)
                i++;
            if (j < knots2.size() - d - 1)
                j++;
        }

        c1->refine_knot(std::move(ins1));
        c2->refine_knot(std::move(ins2));
    }
};

} // namespace alg

} // namespace xi
